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Abstract 

Multivariate global polynomial approximations - such as polynomial chaos 
or stochastic collocation methods - are now in widespread use for sensitivity 
analysis and uncertainty quantification. The pseudospectral variety of these 
methods uses a numerical integration rule to approximate the Fourier-type 
coefficients of a truncated expansion in orthogonal polynomials. For problems 
in more than two or three dimensions, a sparse grid numerical integration rule 
offers accuracy with a smaller node set compared to tensor product approxi- 
mation. However, when using a sparse rule to approximately integrate these 
coefficients, one often finds unacceptable errors in the coefficients associated 
with higher degree polynomials. 

By reexamining Smolyak's algorithm and exploiting the connections be- 
tween interpolation and projection in tensor product spaces, we construct a 
sparse pseudospectral approximation method that accurately reproduces the 
coefficients for basis functions that naturally correspond to the sparse grid in- 
tegration rule. The compelling numerical results show that this is the proper 
way to use sparse grid integration rules for pseudospectral approximation. 
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1. Introduction 

As the power and availability of computers has increased, the profile of 
simulation in scientific and engineering endeavors has risen. Computer simu- 
lations that model complex physical phenomena now regularly aid in decision 
making and design processes. However, the complexity and computational 
cost of the codes often render them impractical for design and uncertainty 
studies, where many runs at different input parameter values are necessary 
to compute statistics of interest. In such cases, designers use a relatively 
small number of high fidelity runs to build cheaper surrogate models, which 
are then used for the studies requiring many model evaluations. 

It is now common to use a multivariate global polynomial of the input 
parameters as the surrogate, particularly when one desires estimates of in- 
tegrated quantities such as mean and variance of simulation results. Addi- 
tionally, the polynomial surrogate is typically much cheaper to evaluate as 
a function of the input parameters, which allows sampling and optimization 
studies at a fraction of the cost. In an uncertainty quantification context - 
where the input parameters often carry the interpretation of random variables 

- this polynomial approximation method appears under the labels polynomial 
chaos [U [2] or stochastic collocation [SJ 0] , amongst others. 

One of the primary disadvantages of the polynomial methods is the rapid 
growth in the work required to compute the approximation as the number 
of model input parameters increases; this generally limits the applicability of 
these methods to models with fewer than ten input parameters. To combat 
this apparent curse of dimensionality, many have proposed to use so-called 
sparse grid methods [5], which deliver comparable accuracy for some prob- 
lems using far fewer function evaluations to build the surrogate. The sparse 
grid is a set of points in the input parameter space that is the union of 
carefully chosen tensor product grids. When the tensor grids are formed 
from univariate point sets with a nesting property, such as the Chebyshev 
points, the number of points in the union of tensor grids is greatly reduced 

- although this nesting feature is not necessary for the construction of the 
sparse grids. The points in the sparse grid can be used as a numerical inte- 
gration rule [6j [7] , where the weights are linear combinations of weights from 
the member tensor grids. Alternatively, the interpolating tensor product La- 
grange polynomials constructed on the member tensor grids can be linearly 
combined in a similar fashion to yield a polynomial surrogate [8], since a 
linear combination of polynomials is itself a polynomial. 
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Another popular polynomial representation employs a multivariate or- 
thogonal polynomial basis. When the coefficients of a series in this basis are 
computed by projecting the unknown function onto each basis, the series is 
a spectral projection or Fourier series [91 [10]; this is also known as the poly- 
nomial chaos expansion in an uncertainty quantification context [TJ [2] . The 
series must be truncated for computation; convergence to the true function 
occurs in the mean-squared sense as one adds more basis polynomials. If the 
integrals in the projections are approximated with a numerical integration 
rule, this method is known as a pseudospectral projection [Hi IT2] . These 
integral approximations only require the simulation outputs evaluated at the 
quadrature points of the input space. 

The question that naturally arises is: Which numerical integration rule 
is appropriate to approximate the Fourier coefficients? Some early attempts 
used Monte Carlo integration [13], but its relative inaccuracies overwhelm 
the spectral accuracy of the truncated Fourier series. Other attempts used 
tensor product Gaussian quadrature rules, but they do not scale to high di- 
mensional parameter spaces due to the exponential increase in the number 
of quadrature points with dimension. The sparse-grid quadrature rules have 
shown promise for retaining the spectral accuracy while alleviating the curse 
of dimensionality. However, in practice this approach produces unacceptable 
errors in the coefficients associated with the higher order basis polynomi- 
als, which forces a much stricter truncation than might be expected for the 
number of function evaluations [T3] . 

This paper presents a sparse pseudospectral approximation method (SPAM) 
for computing the coefficients of the truncated Fourier series with the points 
of the sparse grid integration rule that eliminates the error in the coefficients 
associated with higher degree polynomials. This allows the number of terms 
in the expansion to be consistent with the number of points in the sparse- 
grid integration rule. The key is to separately compute the coefficients of a 
tensor product polynomial expansion for each tensor grid in the sparse grid. 
The linear combination of the tensor weights used to produce the sparse-grid 
integration weights is then used to linearly combine the coefficients of each 
tensor expansion. We show that this method produces a pointwise equiva- 
lent polynomial surrogate to the one constructed from a linear combination 
of tensor product Lagrange polynomials. Therefore error bounds from that 
context can be applied directly. 

Recently, in the context of spectral methods for discretized PDEs, Shen 
and coauthors [15J, [16] proposed and analyzed a closely related sparse spectral 
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approximation using a hierarchical basis of Chebyshev polynomials; the hi- 
erarchical structure results in increased efficiency Their computation of the 
coefficients for the hierarchical basis follow a comparable construction to the 
one we present. However, their focus is on approximating the solution to a 
high-dimensional PDE, as opposed to more general function approximation. 

The remainder of the paper is structured as follows. In Section |2j we re- 
view the relationship between Lagrange polynomial interpolation on a set of 
quadrature points and the pseudospectral approximation for univariate func- 
tions; we then extend this analysis to multivariate tensor product approxi- 
mation. Section [2] closes with a review of Smolyak's algorithm. In Section 
[3j we detail the SPAM for approximating the Fourier coefficients using the 
function evaluations at the sparse-grid integration points followed by some 
interesting analysis results. In Section [4j we present numerical experiments 
from (i) a collection of scalar bivariate functions and (ii) an elliptic PDE 
model with parameterized coefficients. In each experiment, we compare the 
approximate Fourier coefficients from the SPAM with ones computed directly 
with the sparse grid integration rule. Finally we conclude with a summary 
and discussion in Section [Sj 

2. Background and Problem Set-up 

In this section, we briefly review the background necessary to understand 
the SPAM; in particular, we examine the relationship between the Lagrange 
interpolation on a set of Gaussian quadrature points and a pseudospectral 
approximation in a basis of orthonormal polynomials. One purpose of this 
review is to set up the notation, which departs slightly from the notation 
in the disparate references. For the orthogonal polynomials, we follow the 
notation of [TT] . 

Consider a multivariate function / : S — > R, where the domain S C M. d 
has a product structure 

S = Si x • • • x S d . (1) 

Define a d-dimensional point s — (si, . . . , Sj) G S. The domain is equipped 
with a positive, separable weight function w : S — > R + where w(s) = 
Wi(si) • ••Wd{sd) and 

/ s a k w k (s k ) ds k < oo, k = l,...,d, a = 1,2,... (2) 
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The Wk are normalized to integrate to 1, which allows the interpretation of 
w(s) as a probability density function. In general, we consider functions 
which are square-integrable on S, i.e. 

J f{s) 2 w(s)ds < oo. (3) 

Such functions admit a convergent Fourier series in orthonormal basis poly- 
nomials, 

oo oo 

/0) = ^2---^2k,..,i d ni 1 (s 1 )---7i id (s d ) = ^/i7ri(s), (4) 

U = l i<j=l ieN d 

where the equality is in the L 2 sense, i = . . . , i^) is a multi-index, and 

/i = Jj(s)7r i (s)w(s)ds (5) 

is the Fourier coefficient associated with the basis polynomial 7Ti(s). The 
7Ti k (sk) are univariate polynomials in Sk of degree i^ — 1 that are orthonormal 
with respect to u>fc(sfc). In general, a pseudospectral method uses a numerical 
integration rule to approximate a subset of the integrals ([5]); the remaining 
terms are discarded. 

While any square-integrable function admits a convergent Fourier series 
in theory, the polynomial approximation methods perform best on a much 
smaller class of smooth functions; we will restrict our attention to such func- 
tion classes when citing appropriate error bounds. Before diving into the 
multivariate approximation, we first review the univariate case. 

2.1. Gaussian Quadrature, Collocation, Pseudospectral Methods 

Consider the problem set-up above with d — 1. Let 7r(s) = [tti(s), . . . , vr n (s)] T 
be a vector of the first n polynomials that are orthonormal with respect to 
the weight function w(s). The components of 7r(s) satisfy a recurrence rela- 
tionship, which we can write in matrix form as 

sir(s) = Jn(s) + /3 ri+ i7r n+ i(s)e n , (6) 

where e n is an n-vector of zeros with a one in the last entry, and J (known 
as the Jacobi matrix) is the symmetric, tridiagonal matrix containing the 
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recurrence coefficients, 



02 «2 03 



Pn—l &n— 1 Pn 



(7) 



The zeros of vr n+1 (s) generate eigenvalue/eigenvector pairs of J by ([6]), which 
we write 

J = QAQ T , A = diag([Ai,...,A„]), (8) 

where Q(i,j) = 7Tj(A.,-)/||7r(Aj)|| are the elements of the normalized eigen- 
vectors. The zeros Aj of 7r n+ i(s) are the points of the n-point Gaussian 
quadrature rule for w(s); the quadrature weights Uj G M + are given by 



(9) 



k(A;)F 



which are the squares of the first component of the jth eigenvector. A Gaus- 
sian quadrature approximation to the integral is written 



f(s)w(s)ds « W{f) 



3=1 



U 4 



f T v. 



(10) 



The Ug denotes the linear operation of quadrature applied to /; the subscript 
q is for quadrature. This notation will be used later when discussing sparse 
grids. The n-vector f contains the evaluations of f(s) at the quadrature 
points, and the n-vector v contains the weights of the quadrature rule. It 
will be notationally convenient to define the matrices 



P(i,j)=7r<(A 



3)1 



W = diag([ A /z^, . . . , y/hQ) 



(11) 



and note that the orthogonal matrix of eigenvectors Q can be written Q = 
PW. 

The spectral collocation approximation of f(s) constructs a Lagrange 
interpolating polynomial through the Gaussian quadrature points. Since 
the points are distinct, the n — 1 degree interpolating polynomial is unique. 
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We write this approximation Up(f), where the subscript I is for Lagrange 
interpolation, as 

n 

f(s) « Um = £/(A,)400 = f T h». (12) 

i=l 

The parameterized vector l(s) contains the Lagrange cardinal functions 

n — X 

««)= n <«> 

By construction, the collocation polynomial interpolates /(s) at the 

Gaussian quadrature points. 

The pseudospectral approximation of /(s) is constructed by first truncat- 
ing its Fourier series at n terms and approximating each Fourier coefficient 
with a quadrature rule. If we use the n-point Gaussian quadrature, then we 
can write the approximation as 

n 

f( S ) » U^f) = E/iTTiOO = f T 7T( S ), (14) 

1=1 

where fa is the pseudospectral coefficient, 

n 

and the vector f contains all coefficient approximations; the subscript p on 
is for pseudospectral. Note that we have overloaded the notation by defining 
ft as the pseudospectral coefficient ( 15 ), instead of the true Fourier coefficient 
in ([5]). We next state two lemmas about the relationship between the spectral 
collocation and pseudospectral approximations for future reference. 

Lemma 1. The vector of evaluations of f at the quadrature points f is related 
to the pseudospectral coefficients f by 

f = QWf = PW 2 f. (16) 



Proof. This is easily verified by equation (15) using the matrices defined in 
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Lemma 2. The pseudo spectral approximation lAp(f) is equal to the spectral 
collocation approximation Up(f) for all s 6 S. 

Proof. By the uniqueness of the Lagrange polynomial interpolation, we can 
write Pl(s) = Since P = QW" 1 , we have l(s) = WQ T 7r(s). Then 

K(f) = f T i(^) 

= f T WQ T 7r(s) 

= f T 7T( S ) 

as required. □ 

Lemma [2] implies that the pseudospectral approximation I4£(f) inter- 
polates f(s) at the Gaussian quadrature points. However, the equivalence 
expressed in Lemma [2] breaks down in two important cases. When the num- 
ber of terms in the orthogonal series is less than the number of points in the 
quadrature rules, the orthogonal series representation no longer produces the 
same polynomial as the Lagrange interpolant. Also, if a quadrature rule that 
is not the Gaussian quadrature rule is used to approximate the Fourier co- 
efficients, then the discrete Fourier transform from Lemma [T] is no longer 
valid. The latter situation may occur if an alternative quadrature rule holds 
practical advantages over the Gaussian quadrature rule. 

Remark 1. We have restricted our attention to orthonormal polynomials and 
Gaussian quadrature rules for a given weight function. However, transfor- 
mations similar to Lemma^ apply for Chebyshev polynomials and Clenshaw- 
Curtis quadrature rules using a fast Fourier transform. For an insightful 
discussion of the comparisions between these methods of integration and ap- 
proximation, see JT8$ . 

2.2. Tensor Product Extensions 

When d > 1, the above concepts extend directly via a tensor product 
construction. For a given multi- index n = (nj., . . . , rid) £ N d , it is convenient 
to define the set of multi-indices 

X n = : i e N d , 1 < i k < n k , k = 1, . . . , d}. (17) 

We can use this index set to reference components of the tensor product 
approximations. 
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Tensor product Gaussian quadrature rules are constructed by taking cross 
products of univariate Gaussian quadrature rules. For multi-index n, let G n 
be the set of d-variate Gaussian quadrature points, 

G n = {A i = (A n ,...,AJ : iel n }, (18) 

where the points X ik with i k = 1, . . . , are the univariate quadrature points 
for Wk(sk). The associated quadrature weights W n are given by 

W n = {v x = v ix ---v id : i€Z n }. (19) 

In words, the tensor product quadrature weights are products of the univari- 
ate weights. To approximate the integral of f(s), compute 



f(s)w(s) ds^(U^®---®U^)U) 
s 

ni n d 



ii=l ij=l 

where f n is the vector of function evaluations at the tensor grid of quadrature 
points, and u n is a vector of the tensor product quadrature weights. 

The spectral collocation approximation on the points G n uses a basis of 
product-type Lagrange cardinal functions. Define the vector of these basis 
polynomials by 

ln(s) = l ni (si)®---®l nd (s d ), (20) 
where l nk (sk) is a vector of the univariate Lagrange cardinal functions con- 



structed on the univariate quadrature rule defined by X ik ; see (13). Then the 
tensor product spectral collocation approximation for the multi-index n is 
given by 

/(*)«(ar®-"®*C)(/) ( 21 ) 

ni n d 

= J2---J2f(K,--.,K)^(si)---^ d (s d ) (22) 

ii=l i d =l 



(23) 

f„ T ln(*). (24) 
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The tensor product pseudospectral approximation uses a product type mul- 
tivariate orthonormal polynomial basis, which is simply a Kronecker product 
of the univariate orthonormal polynomials. For a multi-index n, let 7r nk {sk) 
be the vector of univariate polynomials that are orthonormal with respect to 
Wfc(sfc) for k — 1, . . . , d. Then the vector 

^n(s) = 7T ni (>i) <g> • • • ® 7T nd (s d ) (25) 

contains polynomials that are orthonormal with respect to w(s); we can 
uniquely reference a component of the vector 7r n (s) by 7Ti(s) with i G X n . 
The tensor product pseudospectral approximation for the multi-index n is 
given by 



{U^®---®U;«){f) (26) 

"1 1l d 

= ■ - ^2 k,...,i d ^h(si) ■ ■ - Tii^Sd) (27) 
iex n 

= (29) 

where f n is the vector of pseudospectral coefficients 

f* = Yl " ' Yl • ■ ■ ' *ii( X h) ■ ■ ■ KuiK) v h ■ ■ ■ v u (3°) 

ii=i jd=i 

= J2f^ n ^»r ( 31 ) 

The extensions of Lemmas [T] and [2] are then straightforward. For the multi- 
index n, define the matrices 

Q = Q ni <8>---<8>Q„ d , P = Pni®-"®Pn d , W = W ni ®---®W nd . 

(32) 

The proofs of Lemmas [T] and [2] hold with 

f = f n , f = f n , l( S ) = l n ( S ), 7T(S) = 7T n ( S ). (33) 

This is easily verified by employing the mixed product property of Kronecker 
products. In words, we have that the Lagrange interpolant constructed on 
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a tensor product of Gaussian quadrature points (i.e., tensor product collo- 
cation) produces the same polynomial approximation as a truncated Fourier 
expansion with a tensor product basis, where the coefficients are computed 
with the tensor product Gaussian quadrature rule (i.e., tensor product pseu- 
dospectral). This equivalence occurs when the number of quadrature points 
in each variable is equal to the number of univariate basis polynomials in 
each variable; in other words, the number of points is the maximum degree 
plus one in each variable. 

2.3. Smolyak's Algorithm and Sparse Grids 

The inescapable challenge for tensor product approximation is the expo- 
nential increase in the work required to compute the approximation as the 
dimension increases. An n-point quadrature rule in each of d dimensions 
uses n d function evaluations. Thus, tensor product approximation quickly 
becomes infeasible beyond a handful of dimensions. Smolyak's algorithm [TH] 
attempts to alleviate this curse of dimensionality while retaining integration 
and interpolation accuracy for certain classes of functions. 

The majority of sparse grid applications in the literature rely on Smolyak's 
algorithm. The most common derivation starts by defining a linear operation 
(e.g., integration, interpolation, or projection) on a univariate function. We 



can generalize the notation used in (10), (12), and (14) by writing the linear 



operation as U m (f). However, it is common to reinterpret the parameter 
m in this context as a choice for how the number of points grows as m is 



incremented. For example, n m = m for m > would correspond to (10) 



(12), and (14). Another common growth relationship is 

n m = 2 m - 1, m > 1. (34) 

Such a relationship is useful when the quadrature /interpolation point sets are 
nested, i.e., the points of the n-point rule are a subset of the points in the 
2n + 1 rule. This notably occurs for rules based on (i) Chebyshev points [TF] , 
(ii) Gauss-Patterson)^] quadrature formulas [20], or (iii) equidistant points. In 
the case of a closed region of interpolation/integration, one may include and 
reuse the endpoints of the interval in the sequence of approximations; see for 



2 The Gauss-Patterson rules contain a specific pattern of nesting that is not applicable 
for arbitrary n. See the reference for further details. 
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example the popular Clenshaw-Curtis integration rules [21J. Nested point 
sets can greatly increase efficiency if f(s) is very expensive. 

Define |m| = mi + - ■ - + md- Given a univariate linear operator, Smolyak's 
method can be written 

A = c ( m ) (^ mi ® " • " ® U md ). (35) 

In the standard formulation [HI |6], the set of admissible multi-indices X is 

X = {meN d : Z + 1 < |m| < Z + d } (36) 
for a given level parameter /. In this case, the coefficients c(m) are 

c (m) ^(-ir-H( ;+ ^ H ). (37) 

However, adaptive and anisotropic versions of Smolyak's algorithm may con- 
tain different choices for X and c(m); such variations are useful if a function's 
variability can be primarily attributed to a subset of the inputs. See [221 E3] 
for details on such methods. 

For our purposes, it is sufficient to note that Smolyak's algorithm amounts 
to a linear combination of tensor product operations. The specific tensor 
products are chosen so that no constituent tensor grid contains too many 
nodes. In the case of nested univariate rules, a node may be common to 
many tensor products. In practice, one may structure the computation to 
evaluate the function once per node in the union of tensor product grids - as 
opposed to once per node per tensor grid. This greatly simplifies the sparse 
grid integration, which can be written as a set of nodes and weights. If a 
node is common to multiple constituent tensor grids, then its corresponding 
weight is computed as a linear combination of the tensor grid weights; the 
coefficients of the linear combination are exactly c(m). It is worth noting 
that the weights of a sparse grid rule can be negative, which precludes its 
use as a positive definite weighted inner product. 

3. Sparse Pseudospectral Approximation Method 

In practice, one may wish to take advantage of the relatively small number 
of points in the sparse grid quadrature rule when computing a pseudospectral 
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approximation. This is often done by first truncating the Fourier series rep- 
resentation of f(s) (see ([5])), and then approximating its spectral coefficients 
with a sparse grid quadrature rule. Unfortunately, choosing the parameters of 
the sparse grid rule that will accurately approximate the integral formulation 
of the Fourier coefficient is not straightforward. The is because - in constrast 
to tensor product approximation - the Lagrange interpolating polynomial is 
not equivalent to a truncated pseudospectal approximation with sparse grid 
integration, where the number of basis polynomials is equal to the number of 
points in the quadrature rule. The general wisdom has been to truncate con- 
servatively for a sparse grid quadrature rule constrained by a computational 
budget; such heuristics become more complicated when anisotropic sparse 
grid rules are used. 

The SPAM approaches this problem from a different perspective; it is 
merely the proper application of Smolyak's algorithm to the tensor product 
pseudospectral projection. We take advantage of the equivalence between 
tensor product pseudospectral and spectral collocation approximations to 
construct spectral approximations that naturally correspond to a given sparse 
grid quadrature rule. In essence, since the sparse grid quadrature rule is 
constructed by taking linear combinations of tensor product quadrature rules, 
we can take the same linear combination of tensor product pseudospectral 
expansions to produce an approximation in a basis of multivariate orthogonal 
polynomials; a linear combination of expansions can be easily computed by 
linearly combining the pseudospectral coefficients corresponding to the same 
basis polynomial. Each tensor product pseudospectral expansion is simply 
a transformation from the Lagrange basis using Lemma [Tj In the numerical 
examples of Section |4| we show compelling evidence that this procedure is 
superior to directly applying the sparse grid quadrature rule to the integral 
formulation of the Fourier coefficients. 

More precisely, let I and c(m) be the admissible index set and coefficient 
function for a given sparse grid quadrature rule. Then the sparse pseudospec- 
tral approximation is given by 



f(s) » A p (f) (38) 

= c ( m ) W 1 ® • • • ® K d )U) (39) 

mel 

= 5>(m)f£7r m ( S ) (40) 
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where f m and 7r m (s) are defined as in (29). In practice, we linearly combine 
the coefficients corresponding to common basis polynomials. With a slight 
abuse of notation, let {tv(s)} be the set of basis polynomials corresponding 
to a vector 7r(s); the common basis set for A p (f) is defined by 

n= (J {*-(«)}■ (4i) 



Then we can write 



Mf)= T, (42) 



w ) - / J 

vr(s)en 



The coefficient corresponding to 7r(s) is given by 



where 



A = ^c(m)/ iim , (43) 



, _ i fi if 7r(s) = 7Ti(s) with i e X m , ^ 
1 '"' ' otherwise. 



In words, (42) simply rearranges the terms in the sum so that each polyno- 
mial basis appears only once. The next theorem allows us to apply existing 
analysis results for sparse grid interpolation schemes to the SPAM. 

Theorem 1. Under the conditions of Lemma \^ the sparse pseudo spectral 
approximation A p (f) is point-wise equivalent to the sparse grid interpolation 
approximation. 

Proof. Using the tensor product version of Lemma [2j we can write 

A(/) = E c ( m ) f - 7r -( s ) 

m£l 

= C ( m )Om(s), 



where f m and l m (s) are defined in (21). This completes the proof. □ 

As a result of this theorem, all of the error analysis for sparse grid colloca- 
tion and interpolation methods applies directly to the sparse pseudospectral 
approximation. We refer the interested reader to references [SI El E] for such 
details. Next, we prove an interesting fact about the mean of A p (f). 
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Corollary 1. The mean of the sparse pseudospectral approximation A p (f) 
is equal to the mean of f(s) approximated with the associated sparse grid 
quadrature rule. 

Proof. By orthogonality, the mean of a polynomial expanded in an orthonor- 
mal basis is equal to the coefficient of the zero degree term, which is 1. Define 
/i to be the coefficient of the constant term in A p (f). The constant term 
also appears in each constituent tensor product pseudospectral approxima- 



tion; denote this by / l m for the multi-index m. Therefore, by (44), 

= J>(m) (ZC®-"®ZC)(/)> 

mel 

which is exactly the definition of sparse grid integration. □ 

3.1. Discrete Orthogonality 

We will see in the numerical results in the next section that - across 
all test cases - the pseudospectral coefficients corresponding to the higher 
order polynomials are inaccurate when computed directly with the sparse 
grid integration rule. This occurs because the higher order basis functions are 
not orthonormal with respect to the sparse grid quadrature rule. However, 
when the integrations are performed using the SPAM, the basis polynomials 
are orthonormal. This becomes apparent by looking at the SPAM coefficients 



for each basis polynomial in the set n from (41). 



Theorem 2. Let f(s) = <p(s) for some 0(s) G II from (41). Then 



U={1 lf l (s) = 0(S) ' (45) 
'0 otherwise, v ' 



where f n is from (43). 



Proof. Using Proposition 3 from [8], we have A p (f) = f for / = e IT, 
which implies that A p is a projector for the polynomial space defined by 
span (n). Noting that the elements of II are linearly independent completes 
the proof. □ 
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(a) SPAM 
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(b) Sparse grid integration 



Figure 1: Orthonormality of the elements in II from (41) using SPAM or 



directly approximated with the sparse grid integration rule with dimension 
d = 2 and level I = 4. The sparse grid was built from univariate Gauss- 



Legendre quadrature rules with growth rule (34). 



Figure la numerically verifies the orthonormality of the elements of II 



using the SPAM; Figure [Lb] demonstrates the loss of orthonormality for the 
higher order elements of II for a discrete inner product defined by the points 
and weights of the sparse grid integration rule. We use d = 2 and I = 4, 
and we order the basis polynomials by their degree. Notice that some of 
the lower order basis polynomials are orthonormal with respect to a discrete 
norm defined by the sparse grid quadrature rule. This is due to the degree 
of exactness of the sparse grid quadrature rule; see [21] for more details. 



4. Numerical Experiments 

In the following numerical experiments, we compare the coefficients com- 
puted with the SPAM to direct approximation of the Fourier coefficients with 
the corresponding sparse grid quadrature rule. To make the comparison fair, 
we apply the sparse grid rule directly to each coefficient corresponding to 
the basis set (41) for the sparse pseudospectral approximation. We con- 



struct each sparse grid rule using (i) univariate non-nested Gauss-Legendre 
quadrature points for a uniform weight function on the square [— 1,1] 2 , (ii) 
n m defined as in (34), and (iii) X and c(m) defined as in (36) and ([37]). 
The choice of the uniform weight function implies the 7Tj(s) are the normal- 
ized Legendre polynomials for the pseudospectral approximation. For all 
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# 


f(si,s 2 ) 


1 

2 




„10„10 


3 
4 


sin(5(si 
1/(2 + 16( 


-0.5)) + cos(3(s 2 - 1)) 
Sl -0.1) 2 + 25(s 2 + 0.1) 2 ) 


5 


(l«i 


-0.2| + |s 2 + 0.2|) 3 



Table 1: The five bivariate test functions. 



experiments, we compute the largest feasible tensor product pseudospectral 
approximation and call the resulting coefficients the truth. In all cases, the 
apparent decay in the tensor product pseudospectral coefficients assures us 
that we have used a sufficiently high order approximation to bestow the title 
truth. 

4-1. Five Bivariate functions 

In the first experiment, we compare both methods on five bivariate func- 
tions; see Table [Tj For each function, we compute a tensor product pseu- 
dospectral approximation of order 255 in each variable - 65,536 total quadra- 
ture points. We plot the log of the magnitude of the pseudospectral coeffi- 
cients with a surface plot to visually observe their decay. We then plot the 
log of the magnitude of the sparse pseudospectral coefficients corresponding 
to a level I = 7 sparse grid compared to the same sparse grid approximation 
of the Fourier coefficients. 

With a level 7 grid, sparse pseudospectral approximation contains a max- 
imum univariate degree of 129. For each test function, the corresponding set 
of figures contains 

(i) the tensor product pseudospectral coefficients up to maximum univari- 
ate degree 100, 

(ii) the SPAM coefficients up to maximum univariate degree 100, 

(iii) the sparse grid integration approximation of the Fourier coefficients up 
to maximum univariate degree 100, 

(iv) a summary plot with coefficients up to univariate degree 129 ordered 
by total order. 

As a general comment, we see that the sparse grid integration produces 
largely incorrect values for coefficients associated with higher degree poly- 
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nomials. More specific comments for the individual test functions are as 
follows: 

1. s\°sl°: This function evaluates the performance of the methods on a 
monomial. We know that coefficients associated with polynomials of 
degree greater than 10 in either si or s 2 should be zero by orthogonality. 
Additionally, since the monomial is an even function over the domain 
with a symmetric weight function, the coefficients corresponding to odd 
degree polynomials in either variable ought to be zero. This is verified 
in the tensor product pseudospectral coefficients and respected by the 
SPAM coefficients. However, the direct sparse integration produces 
non-zero values for coefficients that should be zero. See Figure [2] 

2. e Sl+S2 : This function is analytic in both variables with rapid decay 
of the Fourier coefficients. Observe that the direct sparse integration 
yields large values for coefficients corresponding to the higher order 
polynomials. See Figure [3j 

3. sin(5(si - 0.5)) + cos(3(s 2 - 1)): In the language of the ANOVA de- 
composition |25j, this function has only main effects. Thus, the Fourier 
coefficients for polynomials with mixed terms corresponding to inter- 
action effects should be zero. Again, this is apparent in the tensor 
product pseudospectral coefficients, and it is respected by the SPAM 
coefficients. The direct sparse integration, however, produces non-zero 
values for coefficients of the mixed polynomials; see Figure |4| 

4. 1/(2 + 16(si - 0.1) 2 + 25(s 2 + 0.1) 2 ): The pseudospectral coefficients 
of this rational function decay relatively slowly; notice it needs up to 
degree 40 polynomials in each variable to reach numerical precision, 
according to the tensor product pseudospectral coefficients. The SPAM 
coefficients do a much better job capturing the true decay of the Fourier 
coefficients than the direct integration method, which does not appear 
to decay at all. See Figure |5j 

5. ( | Si — 0.2| + |s 2 +0.2|) 3 : This function has discontinuous first derivatives, 
so we expect only first order algebraic convergence of its Fourier coeffi- 
cients; on a log scale they decay very slowly. However, the interaction 
effects disappear after degree three in either variable. Again, this is 
visible in the tensor product pseudospectral coefficients and respected 
by the SPAM coefficients, but the direct sparse grid integration pro- 
duces non-zero values for coefficients that ought to be zero. See Figure 

El 
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In general, we find that the SPAM coefficients are significantly more accurate 
than the direct application of the sparse grid integration rules to the Fourier 
coefficients. This observation is somewhat counterintuitive. One may expect 
that the sparse grid rule, by evaluating the product of the function times 
the basis polynomial at more locations, would yield a more accurate approx- 
imation. But this is clearly not the case for these examples. The decreased 
accuracy in the coefficients computed with the sparse grid integration rule is 
a result of the nonorthogonality of the basis polynomials with respect to a 
discrete inner product defined by the sparse grid integration rule; see Section 



In Figure 7c, we plot the decay of the truncation error for the sparse 
approximations as the level increases. We approximate the truncation error 
by the sum of squares of the neglected coefficients from the tensor product 
expansion. Since both approximations use the same basis sets, this approxi- 



mate truncation error is identical. In Figures 7a and 7b we plot the decay in 
the error of the approximated coefficients as the level increases; the error in 
the coefficients is squared and summed. We see that the error in the SPAM 
coefficients decays roughly like the truncation error, while the error in the 
direct sparse grid integration does not decay. Of course, this summary plot 
ignores what is most visible in Figures [2]j6j which is that some of the coeffi- 
cients associated with lower order polynomials may be approximated well; it 
is the coefficients of the higher order terms that contain most of the error. 

4-2. PDE with Random Input Data 

The last numerical example we examine is a linear elliptic diffusion equa- 
tion with a stochastic diffusion coefficient. Let D = [0, 1] x [0, 1] and (fl, B, P) 
be a complete probability space. We seek the function u : D x Q — > R such 
that the following holds P-a.e.: 

-V-(a(x,u)VM(x,u)) = l, x E D, 

u(x) = 0, x E dD. ^ ' 

Instead of the whole solution u, we are interested in computing the response 
function 



g(cj) = / u(x,u)dx (47) 
Jd 

which is the spatial mean of u over D. 
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Fi gure 2: Fourier coefficient approximations for sj ^ . 
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Figure 3: Fourier coefficient approximations for e Sl+S2 . 
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Figure 6: Fourier coefficient approximations for (\s\ — 0.2 1 + \s% + 0.2|) c 
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(c) Truncation error 

Figure 7: Comparison of truncation error to error coefficient approximation 
between SPAM and direct sparse grid integration for each of the five test 
functions, numbered according to Table [Tj 
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The diffusion coefficient a(x,u) is modeled as a random field with expo- 
nential correlation: 

C(x, y) = E[a(x, u)a{y, u)] = a 2 e l|x ^ 1|l/L (48) 

where a = 0.1 is the standard deviation of the field and L — 1 is the cor- 
relation length. It is approximated through a truncated Karhunen-Loeve 
expansion 



a(x,u) w a d (x,s(u)) = a (x) + ^ >/\~ k a k (x) s k (u), (49) 



k=l 



where a (x) = \i = 0.2 is the mean of the random field, (A&, a k (x)), k = 
1, . . . , d are eigenvalue-eigenfunction pairs for the covariance operator: 

/ C(x,y)a k (x)dx = X k a k (y), y G D, (50) 

and s = (si, . . . , s<i) are uncorrelated, uniform random variables on [—1, 1]. 
We make the further modeling assumption that the random variables are 
independent. Define T = [— 1, l] d to be the range of s and 

w f 8 ) = { V2 d se{-i,i] d (51) 

\ otherwise 

to be the density of s. The eigenvalues and eigenfunctions are computed using 
a pseudo-analytic procedure described in [2]. The eigenvalues are sorted in 
decreasing order, and we use the first d = 5 eigenvalues/eigenfunctions to 
approximate the random field. 

Let 7Tj : [—1, 1] — > R, % = 1, 2, . . . be the normalized Legendre polynomial 
of order i — 1. For a given multi-index i = (ii, . . . , i^), define the tensor 
product polynomial 

iri{s)=ir il (s 1 )...ir id (s d ). (52) 
Given a set X of multi-indices, we approximate g(u>) by 

9( u ) = ^29iKi(s(uj)) (53) 
iex 

where the unknown coefficients g\ are computed through pseudospectral pro- 
jection using both SPAM and sparse grid integration. For a given s, the 
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corresponding response g is computed by solving 



-V ■ {a d (x, s)Vu(x) = 1, x G D, 
u(x) = 0, x G dD 



(54) 

g = I u{x)dx. 



D 



These equations are discretized using piecewise linear finite elements over 
quadrilateral mesh cells of size 1/512, which gave a spatial error of (9(10 -6 ). 
The resulting linear algebraic equations are solved via preconditioned GM- 
RES using an algebraic multigrid preconditioner with tolerance of 1CT 12 . 
The finite element equations were implemented and solved using a variety of 
packages within the Trilinos solver framework |27j . The resulting SPAM and 
sparse grid integrations were provided by the Dakota package [28] . 



Note that instead of using the growth relationship in (34), we choose 

n m = 2m — 1, m > 1. (55) 
This growth relationship yields tensor grids with many fewer points compared 



to (34). The multiplication factor 2 ensures that all tensor grids will share 
the mid-point of the domain, which reduces the total number of function 
evaluations. The corresponding coefficients of the stochastic response func- 
tion g are plotted in Figure [8] by the degree of the corresponding multivariate 
polynomial. The level parameter for the sparse grid is 4. 

One can see as the order of the polynomials increases, the coefficients 
generated by SPAM decay as they should, whereas those generated through 
direct sparse integration begin to diverge for the higher order polynomials. 
Note, however, that the difference is not as pronounced compared to the 
bivariate test cases. We attribute this to the use of the growth relationship 



(55), as opposed to (34) used with the bivariate functions. 



5. Conclusions 

Sparse grid integration rules are constructed as linear combinations of 
tensor product quadrature rules. By taking advantage of the equivalence be- 
tween the tensor product Lagrange interpolant and a pseudospectral approx- 
imation with a tensor product orthogonal polynomial basis, we can linearly 
combine the tensor product polynomial expansions associated with each ten- 
sor grid quadrature rule to produce a sparse pseudospectral approximation. 
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Figure 8: Comparison of approximate Fourier coefficients of the stochastic 
response (53) of the linear diffusion problem (46) using SPAM and sparse in- 
tegration for dimension d = 5 and level / = 4. The figure plots the coefficients 
according to the degree of associated polynomial. 
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We have numerically compared this approach to direct sparse grid integration 
of the Fourier coefficients. The experiments show convincingly that the di- 
rect integration approach produces inaccurate approximations of the Fourier 
coefficients associated with the higher order polynomial basis functions, while 
the SPAM coefficients are much more accurate. 

The difference between SPAM and the sparse grid integration of the 
Fourier coefficients is present in all Smolyak type approximations - includ- 
ing anisotropic and adaptive variants. While not presented explicitly in this 
paper due to space limitations, the authors have conducted similar studies 
on such variants with similar results. The conclusions are clear. Given a 
function evaluated at the nodes of a sparse grid integration rule, the proper 
way to approximate the Fourier coefficients of an orthogonal expansion is the 
SPAM. 
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